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ABSTRACT 

A microgrid will play an essential role in the future, since it will reduce 
the stress by the increasing electricity demands, as compared to the constraints 
imposed by power plants and transmission lines. Hence, identifying the optimal 
operating point of any distributed generators (DGs) in a microgrid leads them to 
identifying the optimal operation of any generator in therm of cost minimization 
efficiency maximization. 

This thesis proposes a new algorithm for determination of the optimal oper- 
ating point of distributed generators. In this algorithm, the new droop charac- 
teristic slopes (frequency vs. active power and voltage vs. reactive power) of the 
inverter will be detected of the optimized point. Conventionally, inverter-based 
DGs droop characteristic slopes have a constant slope in their maximum and 
minimum capacity in the different operation mode of microgrid (connect to or 
disconnect from the grid). If any change is necessary, the slope will be chosen 
arbitrary. To minimize the deviation of droop characteristic slopes and the other 
objectives (defined by the decision maker), the multi-objective optimization al- 
gorithm can be used. The multi-objective optimization is used in this paper, 



can provide the reasonable decision to determine the new droop characteristic 
slopes in optimized point with minimum deviation. 



This abstract accurately represents the content of the candidate's thesis. I 
recommend its publication. 

Approved: Fernando Mancilla-David 



DEDICATION 

This thesis is dedicated to my wife and son who have supported me all the 
way since the beginning of my studies. They provide me with a great source 
of motivation and inspiration. Also, this thesis is dedicated to my mother and 
memories of my father for instilling in me the importance of hard work and 
higher education. Finally, this thesis is dedicated to all those who believe in the 
richness of learning. 



ACKNOWLEDGMENT 

I would like to express my deepest gratitude to my advisor, Professor Fernando 
Mancilla-David, for his excellent guidance, caring, patience, and providing me 
with an excellent atmosphere for doing research. 

I greatly appreciate Professors Titsa Papantoni, and Professor Dan Connors 
for forming part of my dissertation defense committee. Finally, I would like to 
thank my wife and son; they were always there cheering me up and stood by me 
through the good and bad times. 

I would also like to thank my parents. They were always supporting me and 
encouraging me with their best wishes. 

Finally, I would like to thank my sister in law, Soraya for her supporting. 



CONTENTS 

Figures x 

Tables xi 

Chapter 

1. Introduction 1 

2. Microgrid 3 

2.1 Distributed Energy Resources unit (DERs) 3 

2.2 Microgrid definition and characteristics 3 

2.2.1 Power management strategies for microgrid 6 

3. Optimal Power Flow (OPF) 8 

3.1 Introduction 8 

3.2 Formulation 8 

4. Droop Control 13 

4.1 Concept of droop control 13 

4.1.1 Droop in microgrid 15 

4.2 Load sharing through P-f control 17 

5. OPF with Droop Control 20 

5.1 Mathematic model 20 

5.2 Mathematical formula for optimization 25 

6. Case of Study 30 

6.1 Introduction 30 



vm 



6.2 Microgrid connected to the grid 34 

6.3 Microgrid disconnected from the grid(autonomous mode) 35 

6.3.1 Sharing the power from droop characteristics 35 

6.3.2 Sharing the power with OPF 35 

7. Conclusion 40 

Appendix 

A. "fmincon" in MATLAB 41 

B. The code of software in MATLAB 43 

B.l Reading the file 43 

B.2 Calculate Y-Bus 51 

B.3 Computational code 52 

B.4 Creating fun function 58 

B.5 creatin "noncolon" function 59 

B.6 writing the file after optimization 62 

References 66 



IX 



FIGURES 

Figure 

2.1 Microgrid structure 4 

4.1 Power flowing through a line between two sources 13 

4.2 Frequency droop characteristic 16 

4.3 Voltage droop characteristic 17 

4.4 Frequency droop characteristic for three DGs 18 

5.1 Frequency droop characteristic for DG in bus i 21 

5.2 Voltage droop characteristic for DG in bus i 22 

5.3 Frequency droop characteristics for DG in bus i in different load . 25 

5.4 Voltage droop characteristics for DG in bus i in different load ... 26 

5.5 Direction of frequency droop characteristics changing for DG in bus 

i for different load 26 

5.6 Direction of voltage droop characteristics changing for DG in bus i 

for different load 27 

6.1 The case of study in microgrid 30 

6.2 Case of study with all buses 31 

6.3 Droop characteristic deviation for DG1 38 

6.4 Droop characteristic deviation for DG2 38 

6.5 Droop characteristic deviation for DG3 39 



TABLES 

Table 

6.1 Amount of impedances between busi and bus j in pu 32 

6.2 Data for types of bus 33 

6.3 Data for all of the generators 33 

6.4 Data for all generators and buses in the connected mode with the grid 34 

6.5 Data for microgrid in autonomous mode before optimization .... 36 

6.6 The normalized of objective functions in autonomous mode after op- 
timization 37 



XI 



1. Introduction 

The demand for electricity has increased steadily over the last few years due 
to rapidly population growth. This trend will likely continue over the next few 
decades [18]. Power generation in large centralized power plants such as nuclear, 
fossil fuel (coal, gas), etc. has had an impact on the induced cost. However, the 
customers frequently reside too far from the generating point of electricity. Elec- 
tric power systems and the power grid are currently undergoing restructuring 
due to a number of factors. These factors include: increasing demand, increasing 
uncertainty caused by the integration of intermittent renewable energy sources 
and further deregulation of the industry [5], [9]. 

Due to many economical, practical and technical reasons, the power grid can 
not provide the necessary energy for the consumer. In providing a contingency 
to the power shortage in many areas of demand, it has become necessary to 
search for another way to generate the power required by the consumer. Any 
new decision to produce energy should be researched in parallel with present 
grid supply systems. In the past, all of the power plants have transferred energy 
through transmission lines. However, restrictions in transmission lines force 
the necessities of network design which allows generators to be placed close to 
their loads. In conjunction with the advances in power electronics, the small- 
scaled distributed power generation (non 50-60 HZ) is in the position to change 
the conventional power system [11]. Such change includes the deployment of 
renewable energy sources (RES) such as wind, solar, fuel cells, etc. 



The energy that is generated from natural resources such as wind, sunlight, 
rain, tides and geothermal heat has been named renewable energy (RE), where 
RE unit is an environmentally friendly way of generating power and a serious 
contender to substituting of non renewable energies like coal and oil. Renewable 
energy can be used for a wide range of environments as the resources of such 
environments dictate. Such environments include villages as well as large rural 
areas. Green power refers to a subset of renewable energy that characterized by 
high environmental benefits. 

A large number of research studies have been conducted in many countries 
for substituting conventional energy with RES. For example, It has been re- 
ported that from 1995 to 2000, the U.S. government had spent 1.456 billion 
dollars on renewable energy research [19], while Japan is the rising star of green 
power research. By adjusting their government policies and strong state finan- 
cial support, Japan has made full use of renewable energy. It has also made 
great achievements in wind, solar and ocean energy power generation [20] . Fur- 
thermore, the utilization of green power has disseminated rapidly in Germany 
over the past few years, where it is projected that it will gradually replace the 
nuclear power generation in the next few years. [19]. 



2. Microgrid 

The reduction of fossil fuel sources, pollution of the environment, constraints 
in transmission lines, and the poor efficiency of energy drives us to a new type of 
power generation. These kinds of generators should be located at the distribu- 
tion voltage level and near residential districts. The types of generators in the 
distribution network generate power by using non-conventional energy sources 
(RES). This kind of generator is called a distributed generator (DG) and its 
energy source is called distributed energy resources (DERs). 

2.1 Distributed Energy Resources unit (DERs) 

Distributed generators are small-range power generators that are connected 
to the distribution network, or directly at the customer application [2]. The 
DERs refer to all technologies which can be used to provide energy close to the 
consumer [8]. The high penetration of DGs in the network is able to alleviate the 
pressure on the electric system by feeding into part of the local load. DER units 
contain distributed generator (DG) and distributed storage (DS) with different 
rates of flow [13]. 

2.2 Microgrid definition and characteristics 

In 2001, Professor Robert H.Lasseter from the University of Wisconsin- 
Madison proposed the definition of the microgrid [14] and created a microgrid 
test bed. By this definition, the microgrid has a single unit role in a power 



system. It is defined as the microgrid acting as a group of DER units in the 
same way that distributed voltage can supply electricity and/or heat a group of 
domestic load demand [18]. The proximity of the generator and the load (heating 
and electrical) in the microgrid decrease the capital cost and reduce the losses in 
the transmission lines. A traditional structure of microgrid is shown in Fig. 2.1. 
The microgrid in Fig. 2.1 can operate in two modes: 
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Figure 2.1: Microgrid structure 



1. In this mode, the microgrid is connected to the grid where the grid and 
DER units share the generation of power sufficient for consumers in the 
microgrid section. 



2. In the second mode the microgrid is disconnected from the grid and acts as 



an autonomous island in the microgrid. In this mode all of the DER units 
that exist in the microgrid provide the energy for the consumer. There 
are many reasons for going with this mode. The reasons can range from 
dropped voltage in the grid, a lower price of energy in the microgrid, in 
the event of fault in the grid, and so on. In this case the microgrid can 
stay in stand-alone mode without requiring any power from the grid. 

DER units in microgrid can be connected to the microgrid into two ways: 

1. They are interfaced to the microgrid through rotating machines with re- 
spect to slow dynamics and great inertia. 

2. They are connected through a power electronic converter (PEC) with small 
output impedance and fast dynamics in response to the changing of power. 

So in the transition stage between the two different modes of microgrid 
(stand-alone and connected to the grid), the two connection types can not have 
the same dynamics response. In this situation, in order to increase response time 
for a rotating machine, DS can be used for computational purposes, whereas in 
a steady state condition the rotating machine can provide its power. Discussing 
dynamics in detail is outside the scope of this project. 

DER units and loads in the microgrid can be placed everywhere but not 
necessarily located in the same physical location. All of the type of DER units 
can contribute in the microgrid structure, such as diesel generators, fuel cells, 
and all renewable energy sources, as solar cells (SC), or wind power generators, 
and microturbines as well as other possible types (i.e., geothermal, biomass, 
etc.), as depicted in Fig. 2.1. 



Fig. 2.1 shows the existence of a DS that provides an uninterruptible supply 
of electricity to the microgrid. Currently many types of DS, such as DC bat- 
tery, Flywheel, Supercapacitor (SC), Superconduction Magnetic Energy Storage 
(SMES), etc. are available for the microgrid. 

The microgrid unit has an electrical connection point to the grid that is 
referred to as the point of common coupling (PCC). 

2.2.1 Power management strategies for microgrid 

The objective of this project is to compute the amount of power generated by 
any DER unit in the microgrid. To achieve this objective, a power management 
strategy (PMS) is required. In the microgrid, every DER unit has a different 
rate of flow with no dominant DG in the stand-alone mode of microgrid [12]. 
Therefore, the PMS must consider all types of DER units and create a control 
methodology for energy flow in the microgrid in both modes (connected to and 
disconnected from the grid). In this control, the PMS should supervise the 
active and reactive power flow with high stability and reliability in the system. 

The strategy for control should be to provide the peer-to-peer (P2P) and 
plug-and-play (PnP) concepts in the microgrid. 

In the context of the microgrid, P2P means there is no master control for any 
DER unit, so if one DER unit is lost, the system can continue its operation with 
an extra DER unit. PnP means, any DER units can be added in any location 
of the microgrid without redesigning new control strategies for the system. 

PMS has another responsibility which is to provide the decision in the dy- 
namics in transition between the two modes of the microgrid, but this is outside 



the scope of this project. 



3. Optimal Power Flow (OPF) 

3.1 Introduction 

OPF was presented for the first time by Carpentier in 1962 [4]. OPF is an 

algorithm that optimizes the objective function, such as the total cost of power 
generation, magnitude and phase angle deviation of bus voltage, power losses, 
and emission of pollutants, etc. The objective function can also be predefined 
by the user. OPF finds the optimal set point of the electrical system which 
satisfies the system power flow equations and all of the limit, constraints, and 
contingencies, associated with the power network. 

3.2 Formulation 

In a power grid with N busses, if the voltage magnitude and the phase 

angle of the busses are computed by utilizing the OPF formula then all of the 
variables associated with the system, such as the active and reactive power 
transferred between the buses and also the amount of the current in the lines, 
will be satisfied. In the power system, there are three types of busses: 

1. The PV bus or generator bus, with at least one generator connected to 
the bus. 

2. The PQ bus or load bus, where no generator is connected to the bus. 

3. The slack bus or swing bus, this is defined by the reference for phase angle 
of the voltage and can provide the necessary active and reactive power for 
balancing the system at the optimized solution. 



If we assume the power system includes n busses and n g generators, the 
dimension of the admittance matrix Yb us will be a u n x n" matrix which can be 
denoted as: 
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(3-1) 
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where the elements of y it j, i,j = 1, 2, ...,n of Y& us are constructed as defined in 

[3] as follows: 

• The Yb us is symmetric. 

• The diagonal elements Yu of matrix are equal to the sum of the admittances 
of all the components connected to the ith node. 

• The off-diagonal elements Y^ are equal to the negative of primitive admit- 
tances of all components connected between nodes % and j. 

For this system the complex power injected matrix Sb us is defined as: 



c _ T/ J* 

^bus v bus 1 bus 



(3.2) 



where the Vb us is the vector of the bus voltage and hus is a vector and denotes 
the bus current injected in the bus. Manipulating the equations, the results will 
be as follows: 

With Ifyus = YbusVbus (3-3) 
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(3.4) 


i,j = l,2,...,n 


(3.5) 


i,j = 1,2,.. .,n 


(3.6) 


i,j = 1,2,.. .,n 


(3.7) 


i,j = l,2,...,n 


(3.8) 


i,j = 1,2, ...,n 


(3.9) 



where: 

S'{, MSi : The complex power injection in bus i. 



Sbus z = J~] [\Vi\\Vj\{cosOij + j smOi^jGij - jB tj )] i,j = 1,2, ...,n (3.10) 

n 

Pbus, = ^2 [\Vi\\Vj\(G ij cos 9 i:j + Bij sin 0ij)] i,j = 1,2, ...,n (3.11) 

n 

Qbus, = ^2 [\Vi\\Vj\(G ij sin 9ij -Sy-cos%)] i,j = l,2,...,n (3.12) 

where: 

Pbusi '■ The active power injection in bus i. 
Qbusi '■ The reactive power injection in bus i. 

In the power flow problem, the load is a known variable. In setting up the 
formula, the active and reactive power injection to every bus should consider 
the amount of the load and generation as follows: 

Pbu Si = P 9 i - Pu = f(V, 6) Q^ = Q gi - Q Ll = f(V, 6) (3.13) 
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where: 

P g i : The active power of the generator in the bus i. 
Q g i :The reactive power of the generator in the bus i. 
Pn : The active power of the load in the bus i. 
Qn '■ '■ The reactive power of the load in the bus i. 

According to the optional objective function and its constraints, alternative 
formulas can be presented. The OPF will solve quadratic problem include mini- 
mizing a quadratic function subject to linear and nonlinear equality, inequality, 
that is of interest in this paper. 

The formula in which the objective function minimizes the cost (Ci(P), 



1, 2, ..., n g ) or /and emissions (Tot, 



emmtsstorii j 



1,2, ...,n g ) or/and total losses 



between the lines (Toti osses ) in the system or/and any other objective function 
can be set up as follows: 



Minimize: } ^Cj(P) i = l,2,...,n g 



(3.14) 



i=l 



or/ and 



Minimize: 2_, Tot 



emmi.su/on 



i i = l,2,...,n g 



(3.15) 



i=i 



or/and 



Minimize: Tot 



losses 



(3.16) 



or/and . . . 
subject to: 
inequality constraints: 
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"min S "i S "max 1 — 1, Z, ..., 71 1 ^d.l I J 

Knm < IVJI < V max i = l,2,...,n (3.18) 

*mirii — *gi — *maxi * *-t ^t •••t'lg \0.\.\)j 

|% Fro J < i?ate i,j = 1,2,. ..,n,i ^ j (3.21) 

l-Si^J < iiate i,j = l,2,...,n,i^j (3.22) 

where: 

|SV, From | and |<S!jj T ol are the magnitude of complex power which transferred from 

bus(i) to bus(j) at the sending and receiving point. 

equality constraints: 

n 

P buSi = P gi -P Li = Yl [IMVjKGij cos% + B l3 sin%)] (3.23) 

n 

Qbu Sl = Qgi - Qu = ^2l\Vi\\ Vj | (Gij sin % - Bij cos %)] (3.24) 

3=1 

From equations (3.14) to (3.24), the user can find the variables for the objective 
function. 
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4. Droop Control 

In the power system, active power and frequency are dependent on each 
other. In a rotating machine, increasing the load will result in an increase 
in torque. The torque is related to the rotational speed with relation to the 
frequency of the system. This means that the lower the frequency, the higher 
the load will be. Therefore PMS can try to determine a way to control the 
system with these characteristics when the DER unit in microgrid is interfaced 
to the system with an inverter. 

4.1 Concept of droop control 

To understand the concept of droop control, we return to the power system 
for delivering complex power through a transmission line between two buses. 



Z = R + J XL 



Figure 4.1: Power flowing through a line between two sources 



Fig. 4.1 shows the model of the transmission line between bus % and bus j with 
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impedance Z where: 

Z = R + jXL (4.1) 

and the power flow from the bus % to bus j is defined as: 

Sij = Pij+jQi3 (4.2) 

Sij = VJl 3 (4.3) 

S^V^ (4,) 

by using this definition: 

% = e t - 9j (4.5) 

|t/.|2 _ |y.||y.| e ^j 
S l3 = ^ gl™ (4.6) 

with Euler's formula: 

e 6i ' = cos % + j sin % (4.7) 

we can separate the active and reactive power to this form: 

Pit = X tf + R? ( mi ~ lVjl C ° S% + l^l XLsin %) ( 48 ) 

Qij = xvltf &Wl ~ \Vj\ cos% - \Vj\RBm6a) (4.9) 

If in a transmission line the resistance of the line is much smaller than the 
inductance, then equation (4.8) and equation (4.9) will be changed to this form: 



\ Vi\\Vj \ 
XL 



^ = 4^ sin % ( 41 °) 
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and 

Q ._\^_\}mi coseij (411) 

Furthermore, if the difference in angle 6ij is very small, in computing the 
radian, cos 8ij = 1 and sin#y = 6^. Therefore: 

X LP 

sin% = % = — -^ (4.12) 

I i 1 1 3 I 

and 

\V t \(\V t \-\V,\) = XLQ Z] (4.13) 

From equation (4.12) the difference of the phase angle between two busses (0^) 
is dependent on the active power delivered between the buses, and from equation 
(4.13) the magnitude voltages (|V^| and |V^|) in two buses is influenced by reactive 
power. This implies that by changing 8ij, the active power, and changing the 
\Vi\ and \Vi\, the reactive power can be controlled. 



4.1.1 Droop in microgrid 

In relation to droop control in the microgrid, the frequency is used instead of 
the phase angle for controlling the active power. The microgrid includes two or 
more DGs which use droop characteristics to share the active and reactive power, 
These DGs are interfaced with the inverter. When the load in the microgrid is 
changed, PMS must have specific rules to determine the set point of every DG 
for the sharing of active and reactive power at the optimized point. 

Fig. 4.2 shows the frequency droop characteristic. In this figure the /o is 
the nominal frequency and Pq is the momentary set point for active power of 
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Figure 4.2: Frequency droop characteristic 



the inverter [7]. From the Fig. 4.2, this equation can be writhen as follows: 



fi = fo + KpiP, - P 



(4.14) 



Where K v is the slope of the characteristic and f\ is the the frequency of 
the inverter for generating the power P\ from frequency droop characteristics. 
In Fig.4.3, the voltage droop characteristic is shown. 
From Fig.4.3, the equation can be written as follows: 



\V 1 \ = \V \ + K q (Q 1 -Q c 



(4.15) 



Where Vq is the nominal voltage and Qq is the reactive power of the inverter. 
From equations (4.14) and (4.15) by adjusting P and Q independently, the 
frequency and magnitude of the voltage can be determined [7] . 

4.2 Load sharing through P-f control 
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Figure 4.3: Voltage droop characteristic 

In a microgrid with two or more DGs, changing the load in both modes 
(connectd to the grid or disconnected from grid), PMS uses the frequency droop 
characteristics to change the set-up of the operating point to access the local 
power balance in the new situation of loading. As an example, a microgrid which 
has three DGs with droop characteristics is shown in Fig.4.4. 

Every DG has a different slope in frequency droop characteristics (K pi , K p2 
and K p3 ). 

When the microgrid is connected to the grid, the main grid is large and ro- 
bust enough to compensate any variation in the microgrid [11]. So the frequency 
in the microgrid will be the nominal frequency of the grid (fo). From Fig. 4.4 
and their droop characteristics, the three DGs provide the load Pl as follows: 



Pl 



P,! + P? + P? 



(4.16) 
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Figure 4.4: Frequency droop characteristic for three DGs 



The Pl equals: 



P 



Lv 



±Tj J^nri 



grid 



(4.17) 



The Pgrid denotes the amount of load in the microgrid shared by the grid, 
and the Pl is the total loads in the microgrid. Now if the microgrid for any 
reason is disconnected from the grid, with the same load Pl, the frequency of 
the system will be changed to f\ and the power share for every DG will be P^, 
Pf and Pf, where Pl equals: 



1 i D2 i d3 



In this case: 



p Ll = pi + Pt + p{ 



P Ll = Pl 



(4.18) 



(4.19) 



If the equation (4.14) is developed for every DG, this equation can be written 
as follows: 

fl = f + K Pi (Pi-pi) (4.20) 
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with: 

3 3 

P Ll = J2 P i Pl = J2 P o AP L = P Ll -P Lo (4.21) 

i=l i=l 

Using equation (4.16), equation (4.18), and equation (4.21) the amount of 
power for every DG will be: 

• AP Llh 

P{ = —^ P - + P l (4.22) 
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5. OPF with Droop Control 

An experimental microgrid based on distributed generators {DG) consists 

of several modules: the utility power grid, inverters and loads. This chapter 
will address inverter management and control. Every inverter has autonomous 
droop control characteristics for sharing the active and reactive power amongst 
DGs to meet the load. 

5.1 Mathematic model 

Chapter 4 explained the load sharing process through the p-f controller and 

equation (4.22) shows the amount of power for every DG. In this calculation, the 
total of the power is equal to the loads, but the reality the system is nonlinear 
and the losses should be considered as the equation will be: 



^2 p g t = P Load + Losses (5.1) 



i=i 



The Losses are not constant and are dependent on the voltage of the buses. 

In the autonomous microgrid mode the Losses vary with frequency because the 
impedance is dependent on the frequency. 

To connect many inverters parallel in the microgrid, the frequency of the 
whole system should be by making them the same. The frequency vs. active 
power characteristics for the DG in bus % plot is shown in the Fig. 5.1. 

According to the Fig. 5.1 the following equation for Frequency-droop char- 
acteristic can be written: 
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Figure 5.1: Frequency droop characteristic for DG in bus i 



f = K pi P gi + b. 



(5.2) 



In this equation, the amount of power is dependent on the frequency, so with the 
constant slope in the lower frequency the output power is larger. In stand alone 
mode where the grid is disconnected, every DG should produce more power in 
the same lower frequency. 

The plotting of the voltage vs. reactive power is shown in Fig. 5.2. And 
equation (5.3) shows the relationship of magnitude of the voltage with respect 
to the reactive power. 



w, 



g i 



K qi Q 



qi^Cgi 



'qi 



(5.3) 



Where / and \V g i\ are the frequency and amplitude of the voltage of the DG 
in bus i respectively. K P i and K q i are defined by the corresponding slopes of 
droop. 
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Figure 5.2: Voltage droop characteristic for DG in bus i 

When the microgrid is connected to the grid, the frequency of the system will 
be the same as the utility grid. To share the active and reactive power between 
the DGs at the optimized point, the slopes K pi and K qi can be adjusted. Hence, 
the impedance between two buses % and j can be written as : 



Rij + 27rfLij where: / = 60 



(5.4) 



In autonomous microgrid mode the load will be different, and in order to 
get excellent power sharing, the frequency and magnitude of the DG voltage will 
be changed. So, the change of impedance between two buses i and j caused by 
changing the frequency can be expressed as follows: 

/ 



Xy = 2n60Lij 



and 



Zij(f) — ^ij + JXij 



60 



(5.5) 



For a network with n independent busses, and n g generators, the dimension 
of admittance matrix will be a "n x n" and this matrix will be a function of / 
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as follows (see Chapter 3 and equation 5.5): 



Vn(f) yu(f) --yinif) 

V2l(f) V70.U) ■■V2nU) 



Y bus {f) 



Vnl(f) Vraif) ■ ■ Vnn(f) 



for this network with n buses the following equations can be written: 

q \r t* 

iJbus ' bus 1 bus 

where Ib us is the bus current injection vector. 



b bus — *bus [ *bus \ J ) *bus \ 



^busi 'busi 



^[yijifWbus, 



j=i 



Sbusi = ^2[Vbus t V b * USj y*j(f)] 

3=1 

V; = \VAe j0i V, = \VAe je > 



Y iJ (f) = G lJ (f)+jB tJ (f) % 



, % vj 



Sbus, = J2 [m\VoWHGiAD - JBiAf)]] 

3=1 

where: 

Sbus, '■ The complex power injection in bus i. 



i,j = 1,2,. ..,n 

i,j = 1,2,... ,n 

i,j = 1,2, ...,n 
i,j = 1,2, ...,n 

i,j = 1,2,... ,n 



(5.6) 



(5.7) 

(5.8) 
(5.9) 

(5.10) 

(5.11) 
(5.12) 
(5.13) 



Sbusi = ^[|l / i ||l / j |(cos%+jsin%)(G'i i (/) -jB i:j (f))] i,j = 1,2, ...,n 

(5.14) 
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Pi™ i = ^[\MV j \(G ij (f)co86 ij + B ij (f)sm6 ij )] i,j = 1,2,.. .,n (5.15) 



i=i 



Qtu* =^2[\V i \\V j \(G ij (f)8me ij -B ij {f)cose ij )] i,j = 1,2,... ,n (5.16) 

3=1 



In the power flow problem, the load demands are known variables, so we 
can define the following bus power injections as: 

P buSl = P gi - Pu = f(V, 9, f) Q buSt = Q gi - Q u = f\V, 9, f) (5.17) 

where: 

Pgi : The active power of the DG in the bus i. 
Qgi : The reactive power of the DG in the bus i. 
Pu : The active power of the load in the bus i. 
Qu '■ The reactive power of the load in the bus i. 

From equation (5.2) and equation (5.3) for every DG, the output power (ac- 
tive and reactive) can be computed corresponding to their droop characteristics 

f -b ■ \V -\ — b ■ 

i^ = V Qpi = ^V^ (5 ' 18) 

In optimal power flow the main objective function is the total generation 
cost which is defined as follow: 



Total cost = J2 C *( P ) ( 5 - 19 ) 

where: 
Ci(P) is the cost function for the DG in bus i When the grid is connected to 
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the microgrid, all of the DGs provide their energy froman optimal point at the 
frequency of grid in K pio and K qio . Those are defined by the corresponding slopes 
of their droops. During stand-alone microgrid mode the grid is disconnected 
from the microgrid and the power delivery from each DG will be changed. In 
this variation the amount of power and reactive power sharing will be changed 
according to equation (5.2) and equation (5.3) in the new frequency /. The 
change is shown in Fig. 5.3 and Fig. 5.4. 



Frequency 




p *. P* 



Reactive Power 



Figure 5.3: Frequency droop characteristics for DG in bus i in different load 

However, these points are not optimized points. After calculating the opti- 
mized point, the slope and Y-intersect of the droop will be changed to the new 
value K pi , K qi , b pi and b qi . Fig. 5.5 and Fig. 5.6, depict the direction of changes 
for both droops. 



5.2 Mathematical formula for optimization 

To set up the equations in an optimal power flow, other objectives can be 

added to the optimization problem to minimize the deviation of droop charac- 
teristics. This means the minimizing of AK P i, AK q i,Ab P i and Ab q i. 
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Voltage magnitude 




Reactive Power 



Figure 5.4: Voltage droop characteristics for DG in bus i in different load 



Frequency 



Direction of curve 
can be changed in 
optimization point 




Active Power 



Figure 5.5: Direction of frequency droop characteristics changing for DG in 
bus i for different load 



The equation can be written as follows: 

n 

Minimize: Total cost = J^ d{P) % = 1, 2, ...,n g (5.20) 

and 
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Voltage magnitude 



Direction of curve 
can be changed in 

optimized point 




Reactive Power 



K„.- 



Figure 5.6: Direction of voltage droop characteristics changing for DG in bus 
i for different load 



Minimize: (AK 2 Ab 2 pi , AK 2 Qi , Ab 2 ai 



) i = l,2,...,n g 



(5-21) 



Such that: 



"min — "i _ "max 



v m.ir). 2^ M ^ * rr 



P < P ■ < P 

± mini — 9 % — maxi 



^Cmirii — ^igi _ ^imax. t 

Jmin — J — J'i 



K„ < K„ < K v 

Prntn — l J % — Pmax 



i = 1,2, ...,n 


(5.22) 


i = 1,2, ...,n 


(5.23) 


i = 1,2, ...,n g 


(5.24) 


i = 1,2, ...,n g 


(5.25) 


max 


(5.26) 


i = 1,2, ...,n g 


(5.27) 
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K a . <K a . <K a i = l,2,...,n 

Hrnin — H% — Umax ' ' ' y 

b n . <b r) <b r) i — l,2,...,n 

Prrnn — Pi — Pmax ' ' ' <J 

b qmin < b qi < b qmax i = l,2,...,n g 

\S iJFrom \ < Rate i,j = 1,2, ...,n,i ^ j 
\S iJTo \ < Rate i,j = l,2,...,n,i ^ j 

n 

P buSi = P 9 i -PL i = J2 WWiKGijif) cos % + Bij(f) sin%)] i = 1, 2 

n 

Qbws, = Qgi-QlA = ^2 l\ V i\\ V j\( G ij(f) S ' m0 iJ ~ B ij(f) C0S ^ij)] 1 = 1,2 

3=1 

p / Opi ■ -i r, 

\V -\ —b ■ 

Xqi 
t—^l^-pi -"-pi -^pio 1 ' ' ■■■) ^9 

AK ■ — K ■ — K ■ ?' = 1 2 « 



<-^"pi Dpi "p«o ^ ' ' ■"' ^9 



<—^"qi "qi "qio ^ *-■> ^i •••■> ^g 



5.28) 
5.29) 
5.30) 

5.31) 
5.32) 

...,n 

5.33) 
,...,n 

5.34) 
5.35) 

5.36) 
5.37) 
5.38) 
5.39) 
5.40) 



The objective functions can be evaluated by using the weighted sum method 
that explained in [10], [15], [17], [16], [6]. The weighted objective function for the 
system can be expressed as follows: 

F = W 1 F 1 + W 2 F 2 + ... + W i F i i = l,2,...,4n g + l (5.41) 
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where: F 1 = s ^ j C i (P), 



i=i 



F2 = AK^, F 3 = AK 2 qi , 



F 4 = A^, 


F 5 = Abl 


F 4l _ 2 = AKl, 


F^ = AKl 


F* = Abl, 


Fu+i = Abl 



and, 



W 1 + W 2 + W 3 + ... + W 4i+1 = l (5.42) 
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6. Case of Study 
6.1 Introduction 

The case that we picked up for this study is shown in the Fig. 6.1 [12]. Fij 



Vbase=13.8KV 
Sbase= 100MVA 



Grid 



Utility source 

Q-j) 1000 MWA 

I X/R=22.2 



B-i 



69/13.8 
15MVA 
0.667+j5.33% 



13.8/4.14 
1 .5MVA 



r 

DG3 



? 



L6 



1 .25MVA 
!e + j4B.o% DG2 



13.8/2/4 
3.75MVA 



3 a;:, -is 

.25 MVA 

j.G'j'18.0% 



t 

L3 



13.8/2.4 

2.0MVA 
2.1+J10.94' 



L5 L4 



(^ j 13.8/2.4 

2.5MVA r- p. 

DG1 3 - 2 " ,2 - 3% I 
L2 



-_1_J13.1 



L1 



Figure 6.1: The case of study in microgrid 



6.1 is a single-diagram of a 13.8-KV distribution system used to investigate a 
possible microgrid PMS. This system has three DG units, with different rating, 
ie.,DGl (1.8-MVA), DG2 (2.5-MVA) and DG3 (1.5-MVA). All of the DG units 
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can be dispatched [12]. To study this case, all of the buses for the loads and DGs 
are shown in Fig. 6.2. As seen as in Fig. 6.2, the system has fourteen buses, 



Utility Source 




Figure 6.2: Case of study with all buses 



and three DGs: DG1, DG2 and DG3 connect respectively to Bus6, Busl3, and 
Busl2. The values of the impedance between the lines are stated in Table 6.1. 

This microgrid is connected to the grid through Bus\\. In Table. 6.2 the 
type of the buses with the amounts of the loads (active and reactive power) are 
listed. In the column " Type of bus" of this table, the "No.O" means that there 
is no generator connected to the bus, the "No. 2" means the DG is connected to 
the bus and the "No. 3" means that the bus is connected to the grid. 
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Table 6.1: Amount of impedances between busi and busj in pu. 



Busi 


Busj 


Hij 


XLij 


1 


7 


0.5314 


5.7635 


1 


5 


0.4022 


0.2395 


1 


6 


0.21 


1.094 


1 


14 


0.3976 


0.5127 


2 


14 


0.6141 


0.3066 


3 


14 


0.6065 


1.015 


3 


8 


0.816 


4.8332 


3 


9 


0.6903 


0.000 


14 


11 


0.0818 


0.5826 


12 


4 


0.822 


0.48332 


13 


3 


0.822 


0.48332 


4 


14 


0.3564 


0.2661 


10 


4 


0.822 


4.8332 



The data for all of the generators are listed in the Table 6.3. 
The generator cost curves are presented by quadratic functions [1]. The param- 
eters a, b, and c in Table. 6.3, are the coefficient of the quadratic function and 
are expressed in the following form: 



C t = ai P g 2 t + biP 9i + b t 



1,2, ...,n g 



(6-1) 



Some renewable energy sources such as solar and wind power do not have 
any fuel costs because the resources are free, so in the OPF equation, this type 
of distributed generator resource unit should not be considered for dispatch. For 
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Table 6.2: Data for types of bus 



Bus No. 


Type of bus 


Pl(MW) 


Q L (MVAR) 


Busl 











Bus2 





0.6 


0.42 


Bus3 











Bus4 











Bus5 





0.61 


0.42 


Bus6 


2 








Bus7 





0.600 


0.390 


Bus8 





0.7 


0.45 


Bus9 





0.8 


0.5 


BuslO 





0.9 


0.61 


Busll 


3 








Busl2 


2 








Busl3 


2 








Busl4 


2 









Table 6.3: Data for all of the generators 



DGs.Name 


BusNo. 


Type 


Pmax(Mws) 


Pmin(MW) 


Qmax(MVAR) 


Qmin(MVAR) 


a 


6 


c 


K r 


K 


'<„ 


b. 


DG1 


6 


1 


1.80 


0.000 


1.80 


-1.80 


0.00482 


7.970 


78.00 


-0.0282 


60.0137 


0.0056 


1.0306 


DG2 


13 


1 


2.50 


0.000 


2.50 


-2.50 


0.00156 


7.920 


561.00 


-0.0536 


60.0453 


-0.1488 


1.0183 


DG3 


12 


1 


1.50 


0.000 


1.50 


-1.50 


0.00194 


7.850 


310.00.0 


-0.2537 


60.2417 


-0.2157 


1.0255 


Grid 


11 


1 


1000 


0.000 


1000 


-1000 


0.00139 


7.060 


500.00 


- 


- 


- 


- 



coding in the Table 6.3, if the generator is able to dispatch the type then it is 
(1), otherwise it is (0). Droop characteristic parameters are stated in the last 
four columns of Table 6.3 {K p , b p , K q and b q ). The loads L\ in Bus-?, L 2 in 
Bus 5 , L 3 in Bus2, L± in Bus$, L 5 in Busg and L 6 in Busiq are stated and the 
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value of the loads are expressed in Table 6.2. 

6.2 Microgrid connected to the grid 

In this case, the grid is connected to the microgrid, and the power necessary 
for the loads is shared by DGs and the grid. Table 6.4 shows the values of the 
magnitude, the phase angle of the voltage of each bus, and the amount of the 
power and reactive power shared by any generators, as well as parameters of 
droop characteristics for every DG at the optimized point. 



Table 6.4: Data for all generators and buses in the connected mode with the 
grid 



DGs. Name 


BusNo. 


|VU 


IVte, 


Pl(MW) 


Ql(MVAR) 


P,(MW) 


Q t (MVAR) 


K„ 


b„ 


K, 


b, 


Busl 


- 


1.062 


-0.804 














- 




- 


- 


Bus2 


- 


1.058 


-0.537 


0.600 


0.420 








- 




- 


- 


Bus3 


- 


1.058 


-0.887 














- 




- 


- 


Bus4 


- 


1.063 


-0.585 














- 




- 


- 


Bus5 


- 


1.059 


-0.792 


0.61 


0.420 








- 




- 


- 


Bus6 


DG1 


1.075 


-0.610 








0.5753 


1.1608 


-0.0293 


60.0169 


0.0656 


0.9987 


Bus7 


- 


1.037 


-2.496 


0.600 


0.390 








- 




- 


- 


Bus8 


- 


1.031 


-2.472 


0.700 


0.450 








- 




- 


- 


Bus!) 


- 


1.033 


-2.319 


0.800 


0.500 








- 




- 


- 


BuslO 


- 


1.026 


-2.607 


0.900 


0.610 








- 




- 


- 


Busll 


Grid 


1.066 


0.000 








1.9840 


0.2543 


- 




- 


- 


Busl2 


DG3 


1.073 


-0.660 








0.8708 


0.6927 


-1.0648 


60.9272 


0.3767 


0.8116 


Busl3 


DG2 


1.068 


-1.050 








0.8375 


0.8847 


-0.1078 


60.0903 


0.6452 


0.4973 


Busl4 


- 


1.063 


-0.575 














- 




- 


- 
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6.3 Microgrid disconnected from the grid {autonomous mode) 

If the grid is disconnected from the microgrid, it will not contribute to 

sharing any generated power with the DGs. In this situation, all of the DGs in 

the microgrid should provide the necessary power required for the load without 

the grid. 

6.3.1 Sharing the power from droop characteristics 

Initially, when the grid is disconnected from the microgrid, the droop char- 
acteristics of all inverters will not be changed. New amounts of power generated 
for every DG are computed by equation 4.22 in the same frequency. In Table. 6. 5, 
the amounts of active and reactive power are generated by any DG without any 
change in the parameters of droop characteristic. 



6.3.2 Sharing the power with OPF 

The amount of power for every DG in previous section will not be at op- 
timized point. According to the algorithm expressed in chapter 5, the multi- 
objective function optimization with weighted parameter can be used as follows: 

13 

F = Y,W i F i (6.2) 

i=i 

where: 



F 1 = J2C l (P), 



i=i 
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Table 6.5: Data for microgrid in autonomous mode before optimization 



DGs.Name 


BusNo. 


\v\ m 


O'i,, 


Pl(MW) 


Q L (MVAR) 


P t (MW) 


Q a (MVAR) 


K„ 


K 


K, 


i>. 


Busl 


- 


1.073 


0.125 














- 




- 


- 


Bus2 


- 


1.063 


0.037 


0.600 


0.420 








- 




- 


- 


Bus3 


- 


1.065 


-0.102 














- 




- 


- 


Bus4 


- 


1.067 


-0.007 














- 




- 


- 


Bus5 


- 


1.070 


0.137 


0.61 


0.420 








- 




- 


- 


Bus6 


DG1 


1.091 


1.110 








2.1114 


1.4093 


-0.0293 


60.0169 


0.0656 


0.9987 


Bus7 


- 


1.048 


-1.530 


0.600 


0.390 








- 




- 


- 


Bus8 


- 


1.038 


-1.665 


0.700 


0.450 








- 




- 


- 


Bus!) 


- 


1.041 


-1.514 


0.800 


0.500 








- 




- 


- 


BuslO 


- 


1.031 


-2.009 


0.900 


0.610 








- 




- 


- 


Busll 


Grid 


1.067 


0.000 








0.000 


0.000 


- 




- 


- 


Busl2 


DG3 


1.078 


-0.077 








0.9130 


0.7061 


-1.0648 


60.9272 


0.3767 


0.8116 


Busl3 


DG2 


1.078 


-0.169 








1.2548 


0.9004 


-0.1078 


60.0903 


0.6452 


0.4973 


Busl4 


- 


1.067 


0.000 














- 




- 


- 



2 - 


z ^pv 


4 = 


- ^ 


6 = 


- ^l. 


8 = 


= *& 


10 


= ^ P2 


12 


= A£ 



F 3 

F 7 ~- 
F 9 = 



= *K 



^ii = Ab l 



'13 



*% 



and 



W 1 + W 2 + W 3 + W i + W 5 + W 6 + W 7 + W 8 + W 9 + W 10 + W n + W 12 + W 13 = 1 (6.3) 

For the case study, the problem solved with multi-objective optimization by 
"fmincon" command in MAT LAB (see appendix A). The values of objectives 
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F u F 2 , F 3 , F A) F 5 , F 6 , F 7 , F 8 , F 9 , F 10 , F n , F 12 and F 13 are normalized to be 
consistent with the weights assigned to the objective functions. The functions 
are normalized using the form: 



< 



F - F 



< 1 



(6.4) 



where the Fi is the i th objective function; Fi min is minimum value of i th objective 
function; and Fi max is the maximum value of i th objective function in the different 
weighted, that is assigned to the objective function. Results are shown in Table 
6.6. 



Table 6.6: The normalized of objective junctions in autonomous mode after 
optimization 



Wi 


W 2 :W13 


*i 


F 2 


-Fii 


F 4 


Ff, 


Fa 


F 7 


*8 


F 3 


-Fin 


F n 


F 12 


F 13 


1 





0.085106383 


0.014388834 


0.088864536 


0.697111063 


1.97866E-18 


3.87819E-05 


1.000012322 


1 


0.128438349 


0.124272989 


0.135580785 


1.12556E-17 


1.000139336 


0.9 


0.00833 


0.687943262 


0.032172178 


1.24537E-18 


0.841746065 


0.27036738 


1 


0.096737545 


0.119596111 


0.03213545 


1.000000015 


0.999949561 


0.658812061 


0.788466167 


0.8 


0.0167 


0.709219858 


0.032172178 


0.416655789 


1 


0.142153052 


0.315088821 


0.128846533 


0.581670718 


0.011908363 


0.03949486 


2.45661E-17 


0.781715352 


0.902958985 


0.7 


0.025 


0.680851064 


0.099674757 


0.10919729 


0.480503545 


0.000139537 


3.87819E-05 


0.035064372 


0.418399402 


1.000002212 


0.970045669 


0.078616496 


0.027269824 


0.090394067 


0.6 


0.0333 


0.695035461 


3.92438E-19 


0.164773485 


0.547379823 


1.000019934 


0.215431531 


0.040584768 


-4.67464E-06 


2.70814E-05 


0.049037947 


0.644746978 


0.038824834 


0.210738986 


0.5 


0.04166 


0.758865248 


0.424917685 


1.00003586 


0.196919165 


0.44534645 


0.146438493 


0.013043907 


0.550473308 


-1.50154E-17 


0.013757267 


0.291117184 


0.971427611 


0.166015378 


0.4 


0.05 


0.680851064 


1.000038609 


0.567440913 


0.360510891 


0.051668461 


4.38715E-18 


0.015638986 


0.105922775 


0.009492929 


0.033540092 


0.008405018 


1.000210091 


-1.1165 1E-17 


0.3 


0.0583 


0.737588652 


0.688384534 


0.354241764 


0.218181903 


1.67047901 


0.12638824 


0.011437077 


0.267534592 


0.010167436 


0.003169549 


1.047614681 


0.982890181 


0.000332506 


0.2 


0.0667 


0.723404255 


0.926835785 


0.516317125 


0.118091183 


0.018364031 


0.009272961 


0.003304351 


0.27557615 


0.008642708 


0.01149358 


0.048573027 


0.039976134 


0.018240316 


0.1 


0.075 


0.872340426 


0.37919063 


0.10761944 


0.015645135 


0.100306981 


0.391417765 


1.97157E-06 


0.040762902 


0.000224234 


0.001106104 


0.021129473 


0.014092911 


0.055126289 





0.0833 


1 


0.057976946 


0.018571299 


8.50237E-08 


0.047681697 


0.00042456 


0.005651012 


0.007778609 


0.000230869 


6.6230 1E-07 


0.027824141 


0.008042287 


0.062017075 



The droop characteristics of all DGs are depicted in Fig. 6.3, Fig. 6.4 and Fig. 
6.5 for two conditions: no optimization and optimization with W\ = 0.7, because 
this value of W\ = 0.7 is the most suitable weighted coefficient for optimization. 
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Frequency-droop characteristic for DG1 

W =0.7 

No optimization 



500 1000 1500 2000 2500 3000 
Active power(kw) 
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> 
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0.4 - 
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- - No optimization 
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Figure 6.3: Droop characteristic deviation for DG1 



Frequency-droop characteristic for DG2 



Voltage-droop characteristic tor DG2 



W =0.7 

No optimization 



500 1000 1500 2000 2500 3000 
Active power(kW) 




500 1000 1500 2000 2500 3000 
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Figure 6.4: Droop characteristic deviation for DG2 



From this algorithm, the operator can define the new slope of droop char- 
acteristics at the optimized point. 
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Frequency-droop characteristic for DG3 



Voltage-droop characteristic for DG3 




500 1000 1500 2000 2500 3000 
Active power(kW) 



500 1000 1500 2000 2500 3000 
Reactive power(kVar) 



Figure 6.5: Droop characteristic deviation for DG3 
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7. Conclusion 

This study has presented a methodology for determining optimal power 
flow in a microgrid in a steady state, for either one of two modes of operation, 
i.e., connected to the grid and disconnected from the grid, acting as a stand- 
alone isolated microgrid. This study develops an algorithm which defines the 
droop characteristic parameters of an inverter. The goal is to find the output in 
determining the optimal power to the user. This algorithm assists the user to 
dispatch power in the microgrid. 

If the load changes in the microgrid, this results in changes in the droop 
characteristics in both the slope ( K p and K q ) and the Y-intersect ( b p and 
b q ) deviating at the optimized point. In addition, when the microgrid is dis- 
connected from the grid, the frequency is different from the frequency of the 
grid. From the frequency-droop characteristics the amount of power of DG will 
be changed. The droop characteristic parameters of the inverter-based DGs has 
computed to reach proper load sharing at optimized point rather than to be used 
a fixed value. The power necessary for the system is provided by the optimal 
power flow that minimizes the cost of fuel. 

Transition modes from the connected operating mode to the isolated oper- 
ating mode in a microgrid were not covered in this paper, but only the steady 
state condition in each mode is addressed in this paper. The dynamic transition 
between two modes can be left to a future discussion as an extension of this 
study scope. 
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APPENDIX A. "fmincon" in MATLAB 

fmincon function in MATLAB is proposed to find the minimum of a con- 
strained non-linear multi-variable function. Otherwise fmincon finds a con- 
strained minimum of a function of several variables. The equation can be solved 
by this function expressed as: 
minf(x) 

suchthat : 

c(x) < 

ceq(x) = 

Ax < b 

Aeqx = beq 

lb < x < ub 

x, b, beq, lb and ub are vectors, A and Acq are matrices, c(x) and ceq(x) are 
function that return vectors (©Snonlincon), and f(x) is a function that returns 
a scalar (@fun). f(x), c(x), and ceq(x) can be nonlinear function. The function 
of "fmincon" in MATLAB is expressed as: 

[x, fval, exit flag, output, lambda] = fmincon(@fun, xO, A, b, Aeq, beq, lb, ub, ©nonlcon) 

(Al) 
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where : 
objective 
xO 
A 
b 

Aeq 
beq 

nonlcon 
lb 
up 

X 

fval 
lambda 
exitf lag 
solver 



Objective function 
Initial point for x 

Matrix for linear inequality constraints 
Vector for linear inequality constraints 
Matrix for linear equality constraints 
Vector for linear equality constraints lb Vector 
Nonlinear constraint function 
Lower bound 
Upper bound 

optimized value 

the value of the objective function at the soluti 

Lagrangean multipliers 

specification at algorithm terminatio 

'fmincon ) 
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APPENDIX B. The code of software in MATLAB 
B.l Reading the file 

function S=Read_f ilecf (S) 
f ilename=uigetf ile ( ' * . cf ' ) ; 
S . f ilename=f ilename ; 
f ilecf =fopen(S. filename, 'r') ; 
s=fgetl(f ilecf ) ; 
while strcmp(s, ' ' ) 

s=fgetl(f ilecf ) ; 
end; 

S.CaseName=s; 
while strcmp(s(l:min(3,length(s))) , 'BUS')~=1 

s=fgetl(f ilecf ) ; 
end 

s=fgetl(f ilecf ) ; 
while strcmp(s, ' ' ) 

s=fgetl(f ilecf ) ; 
end; 

S . Bus . number= [] ; 
S.Bus .name=[] ; 
°/ S . Bus . zone= [] ; 
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S.Bus .type=[] ; 
S.Bus.V=[] ; 
S . Bus . ang= [] ; 
S . Bus . MWL= [] ; 
S . Bus . MVARL= [] ; 
S.Bus.MWG=[] ; 
S . Bus . MVARG= [] ; 
S . Bus . Bkv= [] ; 
S.Bus.DV=[]; 
S . Bus . max= [] ; 
S.Bus .min=[] ; 
S.Bus.G=[] ; 
S.Bus.B=[] ; 
while s(l:l)~='-' 

k2=0; 

[kl k2]=find_k(s,k2); 

ite=str2num(s(kl :k2)) ; 
S.Bus. number (ite, l)=ite; 

[kl k2]=find_k(s,k2); 
S . Bus . name=strvcat (S . Bus . name , s (kl : k2) ) ; 

[kl k2]=find_k(s,k2); 
S . Bus . type (ite , 1 ) =str2num (s (kl : k2) ) ; 
[kl k2]=find_k(s,k2); 
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S.Bus.V(ite,l)=str2num(s(kl:k2)) ; 

[kl k2]=find_k(s,k2); 
S . Bus . ang(ite , l)=str2num(s (kl : k2) ) *pi/180 ; 
[kl k2]=find_k(s,k2); 
S.Bus.MWL(ite,l)=str2num(s(kl:k2))/100; 

[kl k2]=find_k(s,k2); 
S . Bus . MVARL ( ite , 1) =str2num (s (kl : k2) ) / 100 ; 

[kl k2]=find_k(s,k2); 
S.Bus.MWG(ite,l)=str2num(s(kl:k2))/100; 
[kl k2]=find_k(s,k2); 
S.Bus.MVARG(ite,l)=str2num(s(kl:k2))/100; 

[kl k2]=find_k(s,k2); 
S.Bus.Bkv(ite,l)=str2num(s(kl:k2)) ; 

[kl k2]=find_k(s,k2); 
S . Bus . DV ( ite , 1 ) =str2num(s (kl : k2) ) ; 

[kl k2]=find_k(s,k2); 
S. Bus. max (ite, l)=str2num(s(kl :k2)) ; 

[kl k2]=find_k(s,k2); 
S.Bus.min(ite, l)=str2num(s(kl :k2)) ; 
[kl k2]=find_k(s,k2); 
S . Bus . G (ite , 1) =str2num(s (kl : k2) ) ; 
[kl k2]=find_k(s,k2); 
S . Bus . B (ite , 1) =str2num(s (kl : k2) ) ; 
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s=fgetl(f ilecf ) ; 
end 



while strcmp(s(l:min(6,length(s))) , 'BRANCH' )~=1 

s=fgetl(f ilecf ) ; 
end 

s=fgetl(f ilecf ) ; 
while strcmp(s,'') 

s=fgetl(f ilecf ) ; 
end; 

S. Branch. fr=[] ; 
S. Branch. to=[] ; 
S. Branch. r=[] 
S. Branch. x=[] 
S. Branch. b=[] 
S . Branch . Mvar= [] ; 
ite=l; 
while s(l:l)~='-' 

k2=0; 

[kl k2]=find_k(s,k2); 

S . Branch . f r ( ite , 1) =str2num (s (kl : k2) ) ; 
[kl k2]=find_k(s,k2); 
S . Branch . to (ite , 1) =str2num(s (kl : k2) ) ; 
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[kl k2]=find_k(s,k2); 
S . Branch . r (ite , 1) =str2num(s (kl : k2) ) ; 
[kl k2]=find_k(s,k2); 
S . Branch . x (ite , 1 ) =str2num (s (kl : k2) ) ; 
[kl k2]=find_k(s,k2); 
S . Branch . b (ite , 1 ) =str2num (s (kl : k2) ) ; 

[kl k2]=find_k(s,k2); 
S . Branch . Mvar (ite , 1) =str2num(s (kl : k2) ) /100 ; 
s=fgetl(f ilecf ) ; 
ite=ite+l; 
end 
while strcmp(s(l:min(9,length(s))) , 'GENERATOR' )~=1 

s=fgetl(f ilecf ) ; 
end 

s=fgetl(f ilecf ) ; 
while strcmp(s,'') 

s=fgetl(f ilecf ) ; 
end; 

if strcmp(s(l) ,'('); 
s=fgetl(f ilecf ) ; 
end; 
while strcmp(s, ' ' ) ; 

s=fgetl(f ilecf ) ; 
end 
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S.G.Bus=[] ; 

S . G . Type= [] 

S.G.Kind=[] 

S . G . Pmax= [] 

S . G . Pmin= [] 

S . G . Qmax= [] 

S . G . Qmin= [] 

S.G.a=[] 

S.G.b=[] 

S.G.c=[] 

'/. S.G.f=[]; 

S.G.Kp=[]; 

S.G.Yp=[]; 

S.G.Kq=[]; 

S.G.Yq=[]; 

S.G.P_ava_UnDis=[] ; 

ite=l; 

while s(l:l)~='-' 

k2=0; 
[kl k2]=find_k(s,k2); 

S . G . Bus (ite , 1) =str2num(s (kl : k2) ) ; 

[kl k2]=find_k(s,k2); 
S . G . Type (ite , 1 ) =str2num (s (kl : k2) ) ; 
[kl k2]=find_k(s,k2); 
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S.G.Kind(ite,l)=str2num(s(kl:k2)) ; 
[kl k2]=find_k(s,k2); 

S . G . Pmax (ite , 1) =str2num(s (kl : k2) ) /100 ; 
[kl k2]=find_k(s,k2); 

S . G . Pmin (ite , 1) =str2num(s (kl : k2) ) /100 ; 
[kl k2]=find_k(s,k2); 

S . G . Qmax (ite , 1) =str2num(s (kl : k2) ) /100 ; 
[kl k2]=find_k(s,k2); 

S . G . Qmin (ite , 1) =str2num(s (kl : k2) ) /100 ; 
[kl k2]=find_k(s,k2); 

S.G.a(ite,l)=str2num(s(kl:k2)) ; 
[kl k2]=find_k(s,k2); 

S.G.b(ite,l)=str2num(s(kl:k2)) ; 
[kl k2]=find_k(s,k2); 

S.G.c(ite,l)=str2num(s(kl:k2)) ; 
°/ [kl k2]=find_k(s,k2); 
% S.G.f (ite,l)=str2num(s(kl:k2)) ; 

[kl k2]=find_k(s,k2); 

S . G . Kp ( ite , 1) =str2num (s (kl : k2) ) ; 

[kl k2]=find_k(s,k2); 

S . G . Yp ( ite , 1) =str2num (s (kl : k2) ) ; 

[kl k2]=find_k(s,k2); 

S.G.Kq(ite,l)=str2num(s(kl:k2)) ; 

[kl k2]=find_k(s,k2); 
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S.G.Yq(ite,l)=str2num(s(kl:k2)) ; 
[kl k2]=find_k(s,k2); 

S . G . P_ava_UnDis (ite , l)=str2num(s (kl : k2) ) /100 ; 
s=fgetl(f ilecf ) ; 
ite=ite+l; 
end 
while strcmp(s(l:min(9,length(s))), , FREQUENCY')~=l 

s=fgetl(f ilecf ) ; 
end 

s=fgetl(f ilecf ) ; 
while strcmp(s, ' ' ) 

s=fgetl(f ilecf ) ; 
end; 
S.f=[]; 
ite=l; 
k2=0; 
[kl k2]=find_k(s,k2); 

S.f=str2num(s(kl:k2)) ; 

S.Bus . pq=f ind(S.Bus.type==0) ; 
S.Bus . pv=f ind(S.Bus.type==2) ; 
S.Bus . s=f ind (S.Bus. type==3) ; 
S.G.Dis=find(S.G.Kind==l) ; 
S.G.UnDis=f ind(S.G.Kind==0) ; 
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'/. S.G.Inv=find(S.G.Type==0); 
N=length(S.G.Bus) ; 
Nl=N-length(S.G.UnDis) ; 
a=max (S . Bus . number) ; 
S.G.Pmax_Dis=S.G.Pmax(find(S.G.Kind==D) ; 

S . G . Pmin_Dis=S . G . Pmin (f ind (S . G . Kind==l ) ) ; 
f close(f ilecf ) ; 



B.2 Calculate Y-Bus 

function S=calculate_YBus(S) 

°/ S.f ilename=' ' ; 

*/. S=Read_f ilecf (S) ; 

S . YB=sparse (diag(S . Bus . G+j *S . Bus . B) ) ; 

S . YB=S . YB+sparse (S . Branch . f r , S . Branch . f r , . . . 
1 . / (S . Branch . r+ j *S . Branch . x) +j *S . Branch . b/2 , 
max (S. Bus .number) , max (S. Bus. number)) ; 

S . YB=S . YB+sparse (S . Branch . to , S . Branch . to , . . . 
1 . / (S . Branch . r+ j *S . Branch . x) +j *S . Branch . b/2 , 
max(S .Bus .number) , max (S .Bus .number) ) ; 

S . YB=S . YB+sparse (S . Branch . f r , S . Branch . to , . . . 
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-1 . /(S. Branch. r+j*S. Branch. x) , max (S. Bus. number) , max (S. Bus. number)) ; 
S . YB=S . YB+sparse (S . Branch . to , S . Branch . f r , . . . 
-1 . /(S. Branch. r+j*S. Branch. x) , max (S. Bus. numb) 



B.3 Computational code 

clear all 
close all 
clc 

°/>starting the reading the data. 
S.f ilename=' ' ; 
S=Read_filecf (S) ; 
°/>For putting the bus that is connected to the grid in the lasr raw: 

S=0rg(S); 
%NN is the number of generators without the grid. 

NN=length(S.G.Bus)-l; 
xl = input ( ' IS THE SYSTEM CONNECTED TO THE GRID ? ' , ' s ' ) 
if (xl=='n') 

S.G.Busl=S.G.Bus(l:NN); 
S.G.Kindl=S.G.Kind(l:NN); 

S.G.Typel=S.G.Type(l:NN) ; 
S . G . Pmaxl=S . G . Pmax ( 1 : NN) 
S . G . Pminl=S . G . Pmin ( 1 : NN) 
S . G . Qmaxl=S . G . Qmax ( 1 : NN) 
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S.G.Qminl=S.G.Qmin(l:NN) ; 

S.G.al=S.G.a(l:NN) 

S.G.bl=S.G.b(l:NN) 

S.G.cl=S.G.c(l:NN) 

'/. S.G.f=[]; 

S.G.Kpl=S.G.Kp(l:NN); 

S.G.Ypl=S.G.Yp(l:NN); 

S.G.Kql=S.G.Kq(l:NN); 

S.G.Yql=S.G.Yq(l:NN); 



S . G . Dis=S . G . Bus (f ind(S . G . Kindl==l) ) ; 

S . G .UnDis=S . G . Bus (f ind(S . G . Kindl==0) ) ; 
S . G . Invl=f ind(S . G . Typel==l) ; 
Nl=length(S.G.Invl); 
°/ Ng is the number of generator, 

°/ Nd is the number of generator that can be dispatchable 
°/ and Nb is the number of busses. 
Ng=length(S.G.Busl) ; 
Nd=length(S.G.Dis) ; 
Nb=max(S. Bus. number) ; 

%S.G.Pmax_Dis is the matrix of maximum reactive power of dispatching 
%generators and S.G.Pmin_Dis is the minimum reactive power of dispatching 
%generators . 
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S . G . Pmax_Dis=S . G . Pmaxl (f ind(S . G . Kindl==l) ) ; 

°/o 

S . G . Pmin_Dis=S . G . Pminl (f ind (S . G . Kindl==l ) ) ; 
S.G.Kp2=S.G.Kpl(S.G.Invl); 
S.G.Yp2=S.G.Ypl(S.G.Invl); 
S.G.Kq2=S.G.Kql(S.G.Invl); 
S.G.Yq2=S.G.Yql(S.G.Invl); 

% Initials value 

xO= [zeros (Nb,l) ;ones(Nb, 1) ; (S.G.Pmin_Dis+S.G.Pmax_Dis) ; (S.G.Qminl+S.G.Qmax 
A=[]; 
b=[]; 
Aeq= [] ; 
beq= [] ; 
°/>upper bond of variable. 

ub= [pi*ones (Nb , 1) ; 1 . l*ones (Nb , 1) ; S . G . Pmax.Dis ; S . G . Qmaxl ; 60 . 8 ; Inf *ones (16,1 

%lower bond of variable . 

lb= [-pi*ones (Nb , 1) ; . 9*ones (Nb , 1) ; S . G . Pmin.Dis ; S . G . Qminl ; 59 . 2 ; -Inf *ones (16 



options=optimset( 'DerivativeCheck' , 'on' , 'Diagnostics' , 'on' , 'Display' , 
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' iter ' , ' FunValCheck ' , ' of f ' , ' GradOb j ' ,'off ,' MaxFunEvals ' , 30000 
'Maxlter ' , 10000 , 'Hessian' , ' of f> ) ; 

[x,fval,exitf lag, output, lambda, grad,hessian]=f mincon (Of un,x0, A, b,Aeq,beq, It 

Nd=length(S.G.Busl)-length(S.G.UnDis) ; 
Ng=Nd+length(S.G.UnDis) ; 
Nb=max(S. Bus. number) ; 
S.Bus.V=x(Nb+l:2*Nb); 
S.Bus.ang=x(l:Nb) ; 

S.Bus.MWG=[] ; 

S . Bus . MVARG= [] ; 

S.Bus.MWG(S.G.Dis)=x(2*Nb+l:2*Nb+Nd)' ; 

S . Bus . MWG (S . G . UnDis) =S . G . P_ava_UnDis (f ind (S . G . Kindl==0) ) ' ; 

S.Bus.MVARG(S.G.Busl)=x(2*Nb+Nd+l:2*Nb+Nd+Ng) ' ; 

f=x(2*Nb+Ng+Nd+l); 

S . G . Kp (S . G . Invl) =x (2*Nb+Ng+Nd+2 : 2*Nb+Nd+Ng+Nl+l) ; 

S . G . Yp (S . G . Invl) =x (2*Nb+Ng+Nd+Nl+2 : 2* (Nb+Nl) +Nd+Ng+1) ; 

S . G . Kq(S . G . Invl) =x (2* (Nb+Nl) +Nd+Ng+2 : 2* (Nb+Nl) +Nd+Ng+N 1+1) ; 

S . G . Yq(S . G . Invl) =x (2* (Nb+Nl) +Nd+Ng+Nl+2 : 2* (Nb+Nl) +Nd+Ng+Nl+N 1+1) ; 

S.f=x(2*Nb+Ng+Nd+l); 
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S . Bus . MWL (f ind (S . Bus . type==2) ) =0 ; 
S . Bus . MVARL (f ind (S . Bus . type==2) ) =0 ; 
else 

S.G.Inv=find(S.G.Type==l) ; 
Ng=length(S.G.Bus); 
Nd=length(S.G.Dis) ; 
Nb=max(S. Bus. number) ; 
S=calculate_YBus(S) ; 

S.G.Dis=S.G.Bus(find(S.G.Kind==l)); 

S . G .UnDis=S . G . Bus (f ind(S . G . Kind==0) ) ; 
NN=length(S.G.Bus)-l; 
S.G.Kpl=S.G.Kp(S.G.Inv); 
S.G.Ypl=S.G.Yp(S.G.Inv); 
S.G.Kql=S.G.Kq(S.G.Inv); 
S.G.Yql=S.G.Yq(S.G.Inv); 

x0= [zeros (Nb,l) ;ones(Nb, 1) ;ones(Nd,l) ; (S.G.Qmin+S.G.Qmax)/2;S.G.Kpl;S.G.Yp 
A=[]; 
b=[]; 
Aeq= [] ; 
beq= [] ; 



ub=[pi*ones(Nb,l) ; 1 . l*ones(Nb, 1) ;S.G.Pmax_Dis;S.G.Qmax; []] ; 
lb=[-pi*ones(Nb,l) ;0.9*ones(Nb, 1) ;S.G.Pmin_Dis;S.G.Qmin; []] ; 
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options=optimset( 'DerivativeCheck' , 'on' , 'Diagnostics ' , 'on' , 'Display' , . . . 
' iter ' , ' FunValCheck' , ' of f ' , ' GradOb j ' , ' of f ' , ' MaxFunEvals ' , 5000 
'Maxlter ' , 10000 , 'Hessian' , ' of f ' ) ; 

[x , f val , exitf lag , output , lambda , grad , hessian] =f mincon (Of unl , xO , A , b , Aeq, beq, 1 



Nl=length(S.G.Inv) ; 

S.Bus.V=x(Nb+l:2*Nb); 

S.Bus.ang=x(l:Nb) ; 

S . G . Kp (S . G . Inv) =x (2*Nb+Ng+Nd+l : 2*Nb+Nd+Ng+Nl) ; 

S . G . Yp (S . G . Inv) =x (2*Nb+Ng+Nd+Nl+l : 2*Nb+Ng+Nd+Nl+Nl) ; 

S . G . Kq(S . G . Inv) =x (2* (Nb+Nl) +Nd+Ng+1 : 2* (Nb+Nl) +Nd+Ng+Nl) ; 

S . G . Yq(S . G . Inv) =x (2* (Nb+Nl) +Nd+Ng+Nl+1 : 2* (Nb+Nl) +Nd+Ng+Nl+Nl) ; 

S . Bus . MWG= [] ; 
S . Bus . MVARG= [] ; 

S.Bus.MWG(S.G.Dis)=x(2*Nb+l:2*Nb+Nd)' ; 

S . Bus . MWG (S . G . UnDis) =S . G . P_ava_UnDis (f ind (S . G . Kind==0) ) ' ; 

S . Bus . MVARG (S . G . Bus) =x (2*Nb+Nd+l : 2*Nb+Nd+Ng) ' ; 

S.f=60; 
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S . Bus . MWL (f ind (S . Bus . type==2) ) =0 ; 
S . Bus . MVARL (f ind (S . Bus . type==2) ) =0 ; 

end 
S=Write_filecf (S); 



B.4 Creating fun function 

function opt=fun2(x,S) 
% opt=0; 

°/. x=ones( 100,1); 
Ng=length(S.G.Busl) ; 
Nd=Ng-length (S . G . UnDis) ; 
Nb=max(S. Bus. number) ; 
Nc=Nb-l; 

Nl=length(S.G.Invl); 
p=100*x(Nc+Nb+l:Nc+Nb+Nd) ; 

kp=x(Nb+Nc+Ng+Nd+2:Nb+Nc+Nd+Ng+Nl+l); 
yp=x(Nc+Nb+Ng+Nd+Nl+2:Nc+Nb+2*Nl+Nd+Ng+l); 

kq=x(Nb+Nc+Nd+Ng+Nl+Nl+2:Nb+Nc+Nd+Ng+Nl+Nl+l+Nl); 
yq=x(Nb+Nc+2*Nl+Nd+Ng+Nl+2:Nb+Nc+2*Nl+Nd+Ng+Nl+Nl+l); 
sl=S.G.a(l:Nd); 
s2=S.G.b(l:Nd); 
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s3=S.G.c(l:Nd); 

c=sl.*(p.~2)+s2.*p+s3; 

Fl=sum(c 
w=l:2*Nl+l; 
w(l)=0.6; 
for i=2:5*Nl+l; 

w(i)=(l-w(l))/(4*Nl); 
end 

% w(2)=l-w(l); 
op=0; 

for ii=l:Nl; 

op=op+w(4*ii-2)*((kp(ii)-S.G.Kp(ii))"2)+w(4*ii-l)*((kq(ii)-S.G.Kq(ii))"2 

end 
opt=op+Fl; 



B.5 creatin "noncolon" function 

function [C, Ceq]=nonlcon(x,S) 
Nd=length(S.G.Busl)-length(S.G.UnDis) ; 
Ng=Nd+length (S . G . UnDis) ; 
Nb=max(S. Bus. number) ; 
Vb=x(Nb+l:2*Nb); 
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ii=l:Nb; 
Vb(ii)=abs(Vb(ii)).*exp(j*x(l:Nb)); 

f=x(2*Nb+Ng+Nd+l); 

Nl=length(S.G.Invl); 

Kp=x(2*Nb+Ng+Nd+2 : 2*Nb+Nd+Ng+Nl+l) ; 

Yp=x (2*Nb+Ng+Nd+Nl+2 : 2*Nb+Ng+Nd+Nl+Nl+l) ; 

Kq=x(2*(Nb+Nl)+Nd+Ng+2:2*(Nb+Nl)+Nd+Ng+Nl+l); 

Yq=x(2*(Nb+Nl)+Nd+Ng+Nl+2:2*(Nb+Nl)+Nd+Ng+Nl+Nl+l); 

% f=x(2*(a+3*N)+l:2*(a+3*N)+N); 

fO=S.f*ones(Nl,l); 

fl=f*ones(Nl,l); 

R=S . Branch . Mvar ; 

SF=Vb (S. Branch. f r).*conj(((Vb(S. Branch. fr)-Vb (S.Branch. to)) ./. . . 

(S . Branch . r+ j *S . Branch . x*f /60) ) +Vb (S . Branch . f r) . * ( j * (S . Branch . b/2) *f /6C 
ST=Vb( S.Branch. to) .*conj ( ((Vb (S.Branch. to) -Vb(S. Branch.fr)) ./. . . 

(S . Branch . r+ j *S . Branch . x*f /60) ) -Vb (S . Branch . to) . * ( j * (S . Branch . b/2) *f /6C 
C=[abs(SF)-R;abs(ST)-R] ; 
C=full(C); 

S . YB=sparse (diag(S . Bus . G+j *S . Bus . B*f /60) ) ; 
S . YB=S . YB+sparse (S . Branch . f r , S . Branch . f r , . . . 
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1 . / (S . Branch . r+ j *S . Branch . x*f /60) + j *S . Branch . b/2*f /60 

max(S .Bus .number) , max (S .Bus .number) ) ; 
S . YB=S . YB+sparse (S . Branch . to , S . Branch . to , . . . 

1 . / (S . Branch . r+ j *S . Branch . x*f /60) + j *S . Branch . b/2*f /60 

max (S. Bus .number) , max (S. Bus. number)) ; 
S . YB=S . YB+sparse (S . Branch . f r , S . Branch . to , . . . 

-1 . /(S. Branch. r+j*S. Branch. x*f/60) , max (S. Bus .number) , max (S. Bus. number)) ; 
S . YB=S . YB+sparse (S . Branch . to , S . Branch . f r , . . . 

-1 . /(S. Branch. r+j*S. Branch. x*f/60) , max (S. Bus .number) , max (S. Bus. number)) ; 
S.YB=full(S.YB); 
Ib=S.YB*Vb; 
Si=Vb.*conj(Ib) ; 
Pi=-S.Bus.MWL; 
Qi=-S.Bus.MVARL; 

Pi (S . G . Dis) =x (2*Nb+l : 2*Nb+Nd) ; 

Pi (S . G . UnDis) =S . G . P_ava_UnDis (f ind (S . G . Kindl==0) ) ; 
Qi(S.G.Busl)=x(2*Nb+Nd+l:2*Nb+Nd+Ng); 
f2=Kp.*Pi(S.G.Busl(find(S.G.Typel==l)))+Yp; 
I Pr=(fl-Yp)./Kp; 

*/. Qr=(abs(Vb(S.Bus.pv))-Yq)./Kq-(S.Bus.V(S.Bus.pv)-S.G.Yq)./S.G.Kq+S.Bus.MV 
% Qr= (abs (Vb (S . G . Bus) ) -Yq) . /Kq; 
V2=Kq.*Qi(S.G.Busl(find(S.G.Typel==l)))+Yq; 
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Ceq= [Pi-real (Si) ; Qi-imag (Si) ; f 2-f 1 ; V2-abs (Vb (S . G . Busl (f ind(S . G . Typel==l) ) ) ) 
% ;Pi(S.G.Bus)-Pr(S.G.Bus);Qi(S.G.Bus)-Qr(S.G.Bus) 
Ceq=full(Ceq) ; 



B.6 writing the file after optimization 

function S=Write_f ilecf (S) 
S . f ilename=uiputf ile ( ' * . cf ' ) ; 
f ilecf =fopen(S. filename, 'w' ) ; 
fprintf (f ilecf , '"/osXn' ,S.CaseName) ; 
fprintf (f ilecf ,' BUS DATA FOLLOWSW); 
for ii=l : length (S. Bus. number) 
fprintf (f ilecf , '°/„4.0f ' ,S. Bus. number (ii)) ; 

fprintf (f ilecf , '°/ s' , S. Bus. name (ii, : )) ; 

°/o 

% fprintf (f ilecf ,'l 1 .Of ' ,S.Bus.area(ii)) ; 
'/. fprintf (f ilecf ,'% 1 .Of ' , S. Bus. zone (ii)) ; 
fprintf (f ilecf , '°/„ l.Of ' , S. Bus .type (ii)) ; 

fprintf (f ilecf ,'% 5.3f ' .S.Bus.V(ii)) ; 

fprintf (f ilecf , "/, 6.3f ' ,S.Bus .ang(ii)*180/pi) ; 

fprintf (f ilecf , '% 7.3f ' , (S.Bus .MWL(ii))*100) ; 
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fprintf (filecf , 7 7.3f ' , (S.Bus .MVARL(ii))*100) ; 

fprintf (filecf , 7 7.4f ' , (S.Bus .MWG(ii))*100) ; 

fprintf (filecf ,7 8.4f , (S.Bus .MVARG(ii))*100) ; 

fprintf (filecf , 7.2 . Of ' ,S . Bus . Bkv(ii) ) ; 

fprintf (filecf , ' %2 . Of ' , S . Bus . DV(ii) ) ; 

fprintf (filecf , 72.0f ' ,S.Bus.max(ii)) ; 

fprintf (filecf , 72.0f ' ,S.Bus.min(ii)) ; 

fprintf (filecf , 7 2. Of ' ,S.Bus.G(ii)) ; 

fprintf (filecf , '% 2. If ' ,S.Bus.B(ii)) ; 

fprintf (filecf , '\n') ; 
end 

fprintf (filecf , ' -999\n' ) ; 
fprintf (filecf /BRANCH DATA F0LL0WS\n'); 
for ii=l : length (S.Branch. r) ; 
fprintf (filecf ,7. 3. Of ' , S.Branch, fr(ii)) ; 
fprintf (filecf , 7„ 3. Of ' , S.Branch. to (ii)) ; 
% fprintf (filecf , 7. 3. Of ', S.Branch. ar(ii)) ; 
% fprintf (filecf , 7 3. Of ', S.Branch. zone (ii)) ; 
% fprintf (filecf , 7 3. Of ', S.Branch. type (ii)) ; 
% fprintf (filecf , 7 3. Of ', S.Branch. circuit (ii)) ; 
fprintf (filecf ,7 5.3f ' , S.Branch. r(ii)) ; 
fprintf (filecf , 7 5.3f ' , S.Branch. x(ii)) ; 
fprintf (filecf ,7 5.3f ' , S.Branch. b(ii)) ; 
fprintf (filecf ,7 5.3f ' , S.Branch. Mvar(ii)*100) ; 
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fprintf (filecf , '\n') ; 

end 

fprintf (filecf , ' -999\n' ) ; 

fprintf (filecf /GENERATOR DATA FOLLOWSW) 

fprintf (filecf , ' (Bus kind Pmax Pmin Qmax Qmin a b c Kp Yp Kq Yq )\n') 

for ii=l:length(S.G.Bus) ; 

fprintf (filecf , 7. 4. Of ' ,S.G.Bus(ii)) ; 
fprintf (filecf . 
fprintf (filecf . 

fprintf (filecf ,7 7.2f ' , (S.G.Pmax(ii))*100) 
7 7.2f ',(S.G.Pmin(ii))*100) 
7 7.2f ',(S.G.Qmax(ii))*100) 
7 7.2f \(S.G.Qmin(ii))*100) 



fprintf (filecf 

fprintf (filecf 

fprintf (filecf 

fprintf (filecf,"/. 7.5f 

fprintf (filecf ,7 6.3f 

fprintf (filecf ,7 6.2f 

fprintf (filecf ,7. 7.4f 

fprintf (filecf ,7. 7.4f 

fprintf (filecf ,7. 7.4f 

fprintf (filecf ,7. 7.4f 

fprintf (filecf ,7 5.2f 

fprintf (filecf , ' \n' ) ; 
end 
fprintf (filecf , ' -99\n ' ) ; 



7. 4. Of ' .S.G.Type(ii)); 
7. 4. Of ' .S.G.Kind(ii)); 



,S.G.a(ii)); 

,S.G.b(ii)); 

,S.G.c(ii)); 

.S.G.Kp(ii)); 

.S.G.Yp(ii)); 

.S.G.Kq(ii)); 

.S.G.Yq(ii)); 

,S.G.P_ava_UnDis(ii)*100) ; 
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fprintf (filecf ,' FREQUENCY IN WHOLE SYSTEMW); 
fprintf (filecf , "/. 7.3f \S.f); 
fprintf (filecf , '\n') ; 
fprintf (filecf , ' -99\n' ) ; 
f close(f ilecf ) ; 
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